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, Abstract 

In this paper we consider an extension of the beta regression model proposed by Ferrari and 
Cribari-Neto (2004). We extend their model in two different ways, first, we let the regression 
D ' structure be nonlinear, second, we allow a regression structure for the precision parameter, moreover, 

this regression structure may also be nonlinear. Generally, the beta regression is useful to situations 
where the response is restricted to the standard unit interval and the regression structure involves 
regressors and unknown parameters. We derive general formulae for second-order biases of the 
maximum likelihood estimators and use them to define bias-corrected estimators. Our formulae 
generalizes the results obtained by Ospina et al. (2006), and are easily implemented by means of 
supplementary weighted linear regressions. We also compare these bias-corrected estimators with 
three different estimators which are also bias- free to the second-order, one analytical and the other 
two based on bootstrap methods. These estimators are compared by simulation. We present an 
■ empirical application. 
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2^ ■ 1 Introduction 
!>■ 
00 

The beta distribution is a very flexible distribution and thus is commonly used to model data restricted 
to some open interval on the Une. The application turns to be more interesting when the interval which is 
being used is the standard unit interval, (0, 1), since the data can be interpreted as rates or proportions. 
' To work with this distribution in a regression manner, several models were defined, see for instance, 

O ! Ferrari and Cribari-Neto (2004), Kieschnick and McCuUough (2003), Paolino (2001), Vasconcellos and 

Cribari-Neto (2005), among others. But the one we will use here is an extension of the one defined 
by Ferrari and Cribari-Neto (2004), mainly because this model is similar to the well known class of 
' generalized linear models (McCuIlagh and Nelder, 1989). Our objective in this paper is to reduce the 

bias of the maximum likelihood estimators (MLEs) for this extended class of beta regression models. 

We extend the model defined by Ferrari and Cribari-Neto (2004) in two ways, namely: (i) We allow 
the regression structure to be nonlinear, such as the exponential family nonlinear models (Cordeiro and 
Paula, 1989) extend the generalized linear models, (ii) We allow a regression structure on the precision 
parameter, this extension is similar to the way the generalized linear models with dispersion covariates 
(e.g., Botter and Cordeiro, 1998) extend the generalized linear models. Further, we allow the regression 
structure on the precision parameter to be nonlinear, an immediate consequence is that we are able to 
model heteroscedasticity in a natural way, by means of this regression structure. Furthermore, we do 
not place any restrictions about the dispersion covariates, in the sense that it can be, and usually is, a 
subset of the covariates of the mean. Moreover, if we see the section on numerical results and the section 
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on application to real data of Ospina et al. (2006), we are able to see that the precision parameter, 
usually has a large variance, thus is not sharply estimated, we then expect the dispersion covariates to 
circumvent this problem. 

The problem of modeling variances has been largely discussed in the statistical literature particularly 

in the econometric area (sec, for instance, Park, 1966; Harvey, 1976). Under normal errors, for instance. 
Cook and Weisberg (1983) and Atkinson (1985) present some graphical methods to detect heteroscedas- 
ticity, Carrol and Ruppert (1988) develop diagnostic procedures using local influence methods for the 
variance parameter estimates in various nonlinear models for the mean, whereas Vcrbyla (1993) com- 
pares full and residual maximum likelihood estimates based on case deletion and likelihood displacement. 
Moving away from normal errors, Smyth (1989) describes a method which allows modeling the dispersion 
parameter in some generalized linear models whereas Bottcr and Cordeiro (1997) present Bartlctt correc- 
tion expressions to improve likelihood ratio test in generalized linear models with dispersion covariates. 
Barroso et al. (2002) present expressions to improve score tests in heteroscedastic Student-t models and 
more recently Taylor and Verbyla (2004) propose a joint modeling of location ans scale parameters in 
Student-t models. Finally, Cysneiros et al. (2007) consider heteroscedastic linear models with symmetric 
errors and discusses the diagnostic aspects of these models. 

It is known that MLEs in nonlinear regression models are generally biased. These bias become 
problematic when the study is being done in small samples. Bias does not pose a serious problem when 
the sample size n is large, since its order is typically 0{n~^), whereas the asymptotic standard error has 
order 0{n^^^^). Several authors have explored bias in regression models. Pike et al. (1979) investigated 
the magnitude of the bias in unconditional estimates from logistic linear models, when the number of 
strata is large. Ratkwosky (1983) uses various examples of normal nonlinear models to relate bias to the 
parameter effects curvature. Also for this class of models, Cook et al. (1986) show that the bias may be 
due to the explanatory variables position in the sample space. Cordeiro and McCuUagh (1991) gave a 
general bias formulae in matrix notation for generalized linear models. 

Bias corrected versions of the MLEs of the parameters that index the beta distribution were obtained 
by Cordeiro et al. (1997) and Cribari-Neto and Vasconcellos (2002). However, such results do not hold 
for models with regression structures. Ospina et al. (2006) obtained expressions for the 0{n~^) bias of 
the parameters of the beta regression models (Ferrari and Cribari-Neto, 2004). 

The method used to obtain expressions of the 0{n~^) bias of the parameters of this class of beta 
regression models is the one given by Cox and Snell (1968). Further, an alternative approach was 
proposed by Firth (1993). He suggested a preventive method of bias reduction by modifying the score 
function prior to obtaining the parameter estimates. It is also possible to perform bias adjustment using 
the estimated bias from a bootstrap resampling scheme, which requires no explicit derivation of the bias 
function. 

The chief goal of this paper is to obtain closed-form expressions for the second order biases of the 
MLEs of the parameters, of the means of the responses, and of the precision parameters of the model. 
The results arc used to define bias corrected estimators to order 0{n~^). We also consider bootstrap 
bias adjustment and the bias adjustment given by Firth (1993). 

The rest of this paper unfolds as follows. In Section 2, we introduce the general class of beta regression 
of interest along with the score function and Fisher's information matrix. In Section 3, we derive a matrix 
expression for the second order biases of the MLEs of the parameters, and consider analytical (corrective 
and preventive) and bootstrap bias correction schemes. We also show how the biases of the MLEs of the 
parameters can be easily computed by means of auxiliary weighted linear regressions. In Section 4, we 
obtain the second order biases of the MLEs of the means of the responses and precision parameters of 
the model. In Section 5, we consider some special cases in detail. In Section 6, we present simulation 
results that show that the proposed estimators have better performance in small samples, in terms of 
bias, than the original MLEs. In Section 7, we consider an empirical example. Finally, the paper is 
concluded in Section 8 with some final remarks. 
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2 The model 



We say that a random variable Y follows a beta distribution with parameters p^q > 0, denoted by 
B{p, q), if the distribution of Y admits the following density with respect to the Lebesgue measure: 

/(y;p,9) = ^J|^2/^-'(i-yr\ ye (0,1), (i) 

where, r(-) is the gamma function. The mean and variance of Y are, respectively 

E{Y) = and Var(r) = ^'^ 



p + q {p + qY{p + q+l)' 

Ferrari and Cribari-Neto (2004) defined a regression structure for beta distributed responses that 
differs from ([T]). Let /x = p/(p + q) and 4> = p + q, i.e., p = ^4> ^nd q = {1 — p)<j3. Under this new 
parameterization, if y ~ B{p, q), then E{Y) — ^ and Var(F) — V{fi)/{1 + 0), where V{iJ,) = /u(l — ^) 
denotes a "variance function", under this parameterization, we will use the notation Y ~ B{^,4>). We 
also note that this parameterization was already known in the statistical literature (see, for instance, 
J0rgensen, 1997, p. 33). Further, (p plays the role of a precision parameter, in the sense that, for fixed 
H, the larger the (/>, the smaller the variance of the response. Using this new parameterization, the beta 
density in ([T]) can be written as 

and the log-density is thus 

log/(y; A*, ^) = \ogT{(b) ~ logr(^(/.) - logr((l - p)^^) + (a*0) logy + {(1 - m)0 - 1} log(l - y), 

with, < /I < 1 and > 0, since p, g > 0. 

Let y — (j/i, . . . , yn)^ be a random sample, where ?/,; ~ B{pi, (f>i), i = 1, . . . ,n. Suppose the mean and 
the precision parameter of yi satisfies the following functional relations: 

giif^t) = Vii = fiixf ; P) and g2{(f>t) =1121^ f2{z'[;0), (3) 

where /3 = (/3i, . . . , /3k)'^ and 9 — {6i, . . . jOh)"^ are vectors of unknown regression parameters which are 
assumed to be functionally independent, (3 € M'^' and 9 S M.^, k + h < n, -qu and r]2i are predictors, and 
ccii, ■ • ■ , Xiq-^ , Zii, . . . , Ziq^ SLie obscrvations on qi and q2 known covariates, which need not to be exclusive. 
We assume that the derivative matrices X = dr]i/d/3 and Z — dri2/d9 have rank k and ft., respectively. 
Moreover, we assume that the link functions gi : (0, 1) ^ M and g2 '■ (0, oo) M are strictly monotonic 
and twice differentiable. A number of different link functions can be used, such as the logit specification 
51 (/i) = log{/x/(l — /z)}, the probit function gi{fJ,) = <i>~^(/i), where $(•) denotes the standard normal 
distribution function, the complementary log-log function gi{p) = log{— log(l — yit)}, among others, and 
for 32, 52(0) = log0, the logarithmic function, (72(0) = the square root function, 52(0) = (with 
special attention on the positivity of the estimates), among others. A rich discussion of link hmctions 
can be found in McCuUagh and Nelder (1989); see also Atkinson (1985, Chapter 7). 
The log-likelihood function for this class of beta regression models has the form 

n 

^(/3,0) =^£,(a.„ (/.,), (4) 

i=l 

where 

= log r((/)i) - log r((i - )(/),,) -f (/Xi 04 - 1) log 

+{{1- fii)^,-l}\og{l~y,); 

lii ~ {flu) T 4'i — g-z^iv^i), as defined in ([3]), are functions of /9 and 9, respectively. It is possible to 
show that this beta regression model is regular, in the sense that all the regularity conditions described 
in Cox and Hinkley (1974, p. 107) hold. It is also possible to show that the MLEs are unique. 
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The components of the score vector, obtained by differentiation of the log-likehhood function with 
respect to the parameters, are given, for r = 1, . . . , /c, as 

where d^ijdrjii = l/g[{^ii), y* = log{yi / {1 - yi)) , fi* = i){iJLi4>i)-'i}j{{l- iJLi)(j)i), and V'(-) is the digammaQ 
function, together with 

Ur{P, 0) = ^ ^ {^^(y* _ + ^(0^) _ ^((1 _ ^^)^^) + log(l _ y,)} 



E ( log ) = - V((l - 



where d(j)i/dri2i — 1/^2 ("^Oi ^'^d i? = 1, . . . , /i. Further, the regularity conditions implies that 

1 - 2/i 
and 

ii;{iog(i - y,)} = ^-((1 - ^^^)^^) - ^-(0.)- 

Consider the complete parameter vector ^ = (/3-^, 6'^)-^. Define the vectors y* = (j/j, . . . , y*)-^, ^* = 
(^1, . . . j^*)"^, w = (ui, . . . ,w„)^, the matrix Ti = diag(d/ii/(i?7H), T2 = diag(d0i/(i?72i), ^ = diag((?!)i), 
with diag(/Xi) denoting the nx n diagonal matrix with typical element = l,...,n, and where 
= Hi{y* — fJ,*) + il^{(j>i) — ~ IJ-i)4>i) + log(l — j/i). Therefore, we can write the (fc + /i) x 1 dimensional 
score vector U{C) in the form (L/^;(/3, 6*)^, t/e(/3, 6')^)'^, with 

U0{f3,9) = X^$Ti(y* -/.*), 
f/e(/3,0) = ~Z^T2v. 

The MLEs of /3 and are obtained as the solution of the nonlinear system ?7(C) = 0. In practice, the 
MLEs can be obtained through a numerical maximization of the log-likelihood function using a nonlinear 
optimization algorithm, e.g., BFGS. For details, see Press et al. (1992). 
Define P as the 2n x (fc + K) dimensional matrix 

P-il »V (0) 



Also, let W be the 2n x 2n matrix 



Z 



W/sg Woe 



with 



Wi3!3 = diag (t)fa. 



W;3e = diag (^<f>i{niai - ip'^l - Hi)(f>i)} 
Wee = diag 



/ d^ij \ I d(jH_ 
\driu) \d7]2i 

, dr]2i J J 

Here, = - Hi)(f>i) + il^' {fii(j>i) and bi = - Hi)(f>i){l - ^i)^ + - Now, using 

(|ni) and ([7]), it is possible to obtain Fisher's information matrix for the parameter vector ( = {/3'^,9^)^ 
as 

K{C) = P'^WP. 



^ We denote generally the polygamma function by ?/)('") (•), m = 0, 1,..., where i/jt™) (x) = (c('"+-'^/dx'"+^) log r(a::), x > 

0. 
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Note that VF^e 7^ 0, thus indicating that the parameters /? and 9 are not orthogonal, in constrast 
to the class of generalized linear models (McCuUagh and Nelder, 1989), where such orthogonality holds. 
Nevertheless, the MLEs C and K{C,) are consistent estimators of C and K{C,), respectively, where K{C,) 
is the Fisher's information matrix evaluated at C,. Assuming that J(C) = lim„^oo K{C,)/n exists and is 
nonsingular, we have that y/n — c}j Nk+hiO, J{C)^^)^ where, denotes convergence in distribution. 
Hence if (r denotes the rth component of it follows that 

(c-c) {Kicrr^/^ ^ N{o,i), 

where K{C,Y'^ is the rth diagonal element of K{C)~'^ . Let K{Q^^ be the {k + i?)th diagonal element of 

K{C)~'^. Then, if < a < 1/2, and represents the 7 quantile of the iV(0, 1) distribution, we have, 

/ - \i/2 . / . \l/2 

for r = 1, . . . ,fc, /3r- ± 91-0/2 (-^(C)'^'^) and On ± qi-a/2 \K{C,)^^\ as the limits of asymptotic 

confidence intervals for (3r and Or: respectively, both with asymptotic coverage of 100(1 — a)%. The 
asymptotic variances of Pr and On are estimated by K{C,Y"'' and K(()^^, respectively. 



3 Bias correction of the MLEs of (3 and 9 

We begin by obtaining an expression for the second order biases of the MLEs of (3 and in this general 
class of beta regression models using Cox and Snell's (1968) general formula. With this expression we 
will be able to obtain bias corrected estimates of the unknown parameters. 

We now introduce the following total log-likelihood derivatives in which we reserve lower-case sub- 
scripts r, s,t,u, . . . to denote components of the /? vector and upper-case subscripts R, S,T^U, . . . for 
components of the vector: U,. = d£/d(3r, Urs = d'^i/ dfirOs, UrsT = d'^£/d(3rd(3sdOT, and so on. 
The standard notation will be adopted for the moments of the log-likelihood derivatives (Lawley, 1956): 
Krs — E(Urs), Kr.s — E(UrUs), Krs,T — ElUrsUx), ctc, where all k's to a total over sample and are, in 

(t) (T) 

general, of order 0{n). We define the derivatives of the moments by Krs = dnrs/dPt, Urs = dnrs/dOr, 
etc. Not all the k's are functionally independent. For example, Krs,t = Krst — f^r) gives the covariance 
between the first derivative of ^(/3, 0) with respect to (3t and the mixed second derivative with respect 
to PnPs- Further, let k^'" = —k^^, k^'^ = —k^'^, k^'^ = —k^^ and k^'^ = —k^^ be typical elements of 
K{Q~^, the inverse of the Fisher's information matrix, which are 0{n^^). 

Let -B(Ca) be the 0{n~^) bias of the MLE for the ath component of the parameter vector ( = 
iCi: ■ ■ ■ Xk,Ck+i, ■ ■ ■ ,Ck+h) = From the general expression for the multiparameter 0{n^^) 

biases of the MLEs given by Cox and Snell (1968), we can write 



B{Ca) 




r,S,U ^ ' B.,S,U ^ ^ 

From ([7]) we note that the entries of the matrix Wf3e are not all zero, which makes the derivation 
cumbersome, since all terms in ^ must be considered. These terms together with the cumulants needed 
to obtain them are given in the Appendix. After some tedious algebra, we arrive at the following 
expression, in matrix form, for the second order bias of /?: 

B0) = KP^X^[KhPpfi + {M2 + Mz)Pf3B + MM + K^''X^[N2G-NiF] 
+ rP'Z^ [M2P/313 + (M4 + M5)Pi30 + MePee] + K^'" [N2F - N^G] , 
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where Mi to AIq and Ni to N3 are given in (HOI)- (123) in the Appendix, Kf^f^, K^^ and K'^'^ are the 
matrices formed by the (r, s)th, (r, S')th and {R, S)th. elements of the inverse of the Fisher's information 
matrix with r,s — 1, . . . , fc, R, S = k + 1^ . . . , k + h, respectively, -P/3/3, Ppe and Pgg are the n x 1 
dimensional vectors containing the diagonal elements of XK^^X'^ , XK^^Z^ and ZK^^Z^ , respectively, 
F — diag(Fi, . . . , Fi = ti-{XiK^^), Xi is a, k x k matrix with elements ff^rju/ dfirl^s, 1 is an n x 1 

vector of ones, G = diag(Gi, . . . , G„)l, Gi = tr{ZiK^^), Zi is a hxh matrix with elements d'^r]2i/d6iid6s. 
We define the 2n x 1-vectors ivi and uj2 as 

^ . ' + (Af2 + M3)Pf}e + M5P00 

' M2Pf3p + {Mi + M5)Ppe + MeP99 



and 

UJ2 



N2G - NiF 
N2F - N3G 



We also consider the k x {k + h) upper block of the matrix K{() ^ given by 

The 0(n~^) bias of /3 can now be written as 

bCP)^K^*P'^{uJi+u^2). (9) 

We now turn to the 0(n^^) bias of 9. Analogously, we have the following expression, in matrix form, 
for the second order bias of 6: 

B{e) = Rf^'^X^ [MiPpp + {M2 + Kh)Ppg + NhPgg] + RI^^ X'^ [N2G - N^F] 
+ i^^^Z^ [M2PPP + (Af4 + Af5)P/3e + MM + K'">Z'^ [N2F - A^gG] , 

Then, considering the h x {k + h) lower block of the matrix R{Q^^, given by 
we can write the 0{n~^) bias of 6 as 

B{e)^R''*P^{uji+U2). (10) 

Thus, combining ^ and pUj) . we are able to find that the 0{n^^) bias of the MLE of the joint vector 
C = (/3^, 0"^)^ can be written, in matrix form, as 

B(C) = R{0-^P'^{uji + UJ2) = [P^WP)-^P'^{uJi + UJ2). 
Now, let — W~^uJi and ^2 ~ W^^lu2, then, the previous expression becomes 

B{c)^{p^wpy^p^w{^i+^2). (11) 

Thus, the 0{n~'-^) bias of C (fTT|) is easily obtained as the vector of regression coefficients in the formal 
weighted linear regression of ^ = ^1 + ^2 on the columns of P with W as weight matrix. 

The ©(n"^) bias pT|) is expressed as the sum of two quantities: (i) Bi = (^P^WP) P^W^i, 
the bias for a linear beta regression model with dispersion covariates with model matrices X and Z, 
and thus generalizes the expression obtained by Ospina et al. (2006), and (ii) an additional quantity 
B2 — [P'^WP^ P^W^2 due to the nonlinearity of the functions /i(xi;/3) and f2{zi;9), and which 
vanishes if both /i and /2 are linear in /3 and 9, respectively. 

Now we can define our first bias-corrected estimator C as 

C = C - B{C), 

where B{Q denotes the MLE of B{(), that is, the unknown parameters are replaced by their MLEs. Since 
the bias B{Q is of order 0{n~'^), it is not difficult to show that the asymptotic normality ( C ^ C ) ^ 



6 



Nk+h{0,J ^(C)) still holds, where, as before, we assume that J(C) — lim„^oo -f^^(C)/" exists and is 

( . 1/2 

nonsingular. From the asymptotic normality of we have that ± 91-0/2 \ K{(^)"-'^ > , for a = 

1, . . . ,k, k + 1, . . . , k + h. The asymptotic variance of is estimated by i4r(C)'"', where if (C)'^" is the ath 
diagonal element of the inverse of the Fisher's information matrix evaluated at 

We now turn to the bias-correction approach developed by Firth (1993) called the "preventive" 
method. This method also remove the 0(n~^) bias and consists of modifying the original score function: 

U*{C)^U{C)+K{C)B{0, 

where K{t^) is the information matrix and i?(C) is the 0{n^^) bias. The solution ^ of J7* = is a 
bias-corrected estimator, to order 0{n~^), for (. For our general class of beta regression models, the 
substitution of B{() by (fTTj) yields the following form for the modified score function: 

U*{C) = U{C)+P^iuJi+L02). (12) 

The estimator obtained as the root of the modified score function in ([T^ . is consistent and 

asymptotically normal: y/n {C ~ C) ~^ -^fc+/i(0, J^^(C)), with J(^) as given before. We also have that 
1/2 

(a ± qi~a/2 {^(C)""} I for a = 1, . . . ,k,k + 1, . . . ,k + h. The asymptotic variance of (a is estimated 
by K{()"-°', where K{()"-'^ is the ath diagonal element of the inverse of the Fisher's information matrix 
evaluated at C- 

A third, and the last approach we consider here, to bias-correcting MLEs of the regression parameters 
is based upon the numerical estimation of the bias through the bootstrap resampling scheme introduced 
by Efron (1979). Let y = (yi, . . . ,yn)^ be a random sample of size n, where each element is a random 
draw from the random variable Y which has the distribution function F = F{Q. Here, C is the parameter 
that indexes the distribution, and is viewed as a functional of F, i.e., C — t{F)- Finally, let C be an 
estimator of ( based on y; we write C = •s(y). 

The application of the bootstrap method consists in obtaining, from the original sample y, a large 
number of pseudo-samples y* = (y^, . . . , y*)^, and then extracting information from these samples 
to improve inference. Bootstrap methods can be classified into two classes, depending on how the 
sampling is performed: parametric and nonpar ametric. In the parametric version, the bootstrap samples 
are obtained from F((), which we shall denote as F^, whereas in the nonparametric version they are 

obtained from the empirical distribution function F, through sampling with replacement. Note that the 
nonparametric bootstrap does not entail parametric assumptions. 
Let Bp{(, C) be the bias of the estimator C, = s{y), that is, 

Bp{C, = Ej.[C - C] = Ep[s{y)] - t{F), 

where the subscript F indicates that expectation is taken with respect to F. The bootstrap estimators 
of the bias in the parametric and nonparametric versions are obtained by replacing the true distribution 
F, which generated the original sample, with F^ and F, respectively, in the above expression. Therefore, 
the parametric and nonparametric estimates of the bias are given, respectively, by 

BF^{C,O^EF^Hy)]-t{F^) and i?^(C, C) = E^[s(y)] - t(/'). 

If B bootstrap samples (y*'^, y*^, . . . , y*^) are generated independently from the original sample y, and 
the respective boostrap replications (C*^, C*^) • ■ • ) C*'^) are calculated, where C*** = s{y*^), b = 1,2, . . . ,B, 
then it is possible to approximate the bootstrap expectations Ef^. [s(y)] and E^[s(y)] by the mean 

^*( ) — ^ ^b=i C***- Therefore, the bootstrap bias estimates based on B replications of ( are 

Bp^iLO^C^-^ ~s{y) and B^(C, C) = C^'^ - s(y), (13) 

for the parametric and nonparametric versions, respectively. 

By using the two bootstrap bias estimates presented above, we arrive at the following two bias- 
corrected, to order 0{n~^), estimators: 

Ci = s(y)-Bf^-(C,C)-2C-C*^-\ 
C2 = s(y)-B^(C,C)-2C-r^-^. 
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The corrected estimates C,i and were called constant-bias-correcting (CBC) estimates by MacKinnon 
and Smith (1998). 

Since we are dealing with regression models and not with a random sample we need some minor 
modifications to the algorithm given above. 

For the nonparametric case, assume we want to fit a regression model with response variable y and pre- 
dictors xi, . . . , Xq^ , zi, . . . , Zq^. We have a sample of n observations pf = {yi,Xii, . . . , Xiq^ , z^i, . . . , ztq^)^ 
i = 1, . . . ,n. Thus we use the nonparametric bootstrap method described above to obtain B bootstrap 
samples of the pf , fit the model and save the coefficients from each bootstrap sample. We can then 
obtain bias corrected estimates for the regression coefficients using the methods described above. This 
is the so-called Random- a; resampling. 

For the parametric case, assume we have the same model as for the nonparametric case, we thus 
obtain the estimates fii and (jji (such as in our case where the distribution is indexed by /i and (jj) 
and using the parametric method described above, we obtain B bootstrap samples for iji from the 
distribution F{jii,(j)i), i — 1, ... ,71. We would then regress each set of bootstrapped values yl on the 
covariates xi, . . . , Xq-^ , zi, . . . , z^^ to obtain bootstrap replications of the regression coefficients. We can, 
again, obtain bias corrected estimates for the regression coefficients using the methods described above. 
This method is called Fixed-x resampling. 



4 Bias correction of the MLEs of /i and 

In this Section we obtain the results that are the most valuable to the practioners, namely, the 0{n~^) 
bias of /X and of (p, since, for practioners, the interest in a data analysis relies on sharp estimates of 
the responses and of the precision parameters. The fact that these results must be computed apart 
comes from the fact that if /3 and 9 are bias-free estimators, to order 0{n~^), it is not true, in general, 
that Hi = gi^{fi{xi;f3)) and 4>i = g2^{f2{zi;9)) will also be bias-free to order 0{n~^). Nevertheless, for 
practioners, it is even more important to correct the means of the responses and the precision parameters 
than correcting the regression parameters. Moreover, the 0{n~^) bias of jj. for the linear beta regression 
model was not presented in Ospina et al. (2006). 

We shall first obtain the 0{n~^) bias of the MLEs of 771 and 772- Using ^ we find, by Taylor 
expansion, that to order 0{n^^): 

h{xj- 0) - h{xf- p) = ^p{vuf0 + i(/3 - PfMP - /3), 

and 

/2(^f ; 9) - /2(zf ; 9) = Veim^fi9 ~ 9) + ^{9 ~ 9f Z,{9 - 9), 

where V^(r/ii) is a fc x 1 vector with the derivatives drju/dPr, ^e(jl2i) is a ft. x 1 vector with the derivatives 
dT]2i/d9ji, and all the other quantities were previously defined in page [HI 

Thus, taking expectations on both sides of the above expression yields to this order 

B(r;i) = XB0) + If, 

and 

B{yn) = ZB{9) + Ig, 

where, F and G were defined in page[Sl and since K^^^ and K^^ are the asymptotic covariance matrices 
of (3 and 9, respectively. 

From similar calculations we obtain to order 0{n^^) 

and 

B{^,)^B{f,2^)^ + \w^r{m^)^. 

ami 2 dr]^i 
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Let Ti and T2 be as in pagejH further, let 5*1 — diag{d'^ fii / dr]f^) and 5*2 — 6ia,g{dP (jii / drjli) ■ Then, we 
can write the above expressions in matrix notation as 

B{fi) = \t^{2XB0) + F) + \s,Ppp (14) 

and 

B{^) = ^T2i2ZB{9) + G) + ^S2Pee, (15) 

where Ppf} and P0O were defined in page [6] and since the asymptotic covariance matrices of r/i and 772 are 
XRl^l^X^ and ZK^^'^Z'^, respectively. 

If we combine (|14p and (fT5| with ^ and (fTO|) , we will have the following explicit expressions for the 
0{n^^) biases of jl and 0, respectively: 

Bi(A) = ^m2XK^*P^iLOi+L02)+F) + ^SiPf3^ 

and 

bS) = iT2(2Zi^'^*F^(wi + CJ2) + G) + ^S2Pee- 

Further, let $ and be the estimators obtained as the root of the modified score function (fT^ . Then 
we also obtain the following formulae for the 0{n~^) biases of fi and (j), respectively: 

B2(fi) = It,{2X0 -P)+F) + ^S.Ppp 

and 

B2i$) - lT2i2Z{§ -d) + G) + y2Pe9. 

Lastly, we can use the bootstrap-based 0{nr^) biases to define, bias corrected estimators of fi and 
<j) to this order. Then, let Bp-{P) be the vector formed by the first k elements of the vector Bp-{t^,(^) 

defined in equation p^ . Bp- (9) be the vector formed by the last h elements of the vector Bp- {(, (), and 
define Bp[(3) and Bp{6) analogously from the vector Bp{C, C) also in equation Thus, we have the 
following alternative expressions for the 0{n^^) biases of // and 0, respectively: 

B^iP) = \t^{2XBp.0) +F) + \SiP0i3 and B^ifi) - ]^T^{2XBpCp) + F) + ]^SiPpp, 

and 

bS) = \t2{2ZBp-{6) + G) + ]^S2Pee and bS) = \T2{2ZBp{6) + G) + ]^S2Pee. 
Therefore, we are now able to define the following second-order bias-corrected estimators for {1 and 

fi^ fi- Biifl), fi^ js - B2{jl), ^ fl - Bsip,) and ^2=A~^4(A) 

and 

4> = $- bS), = 0- M^), ^1=0- bS) and 02 = - ^4(0), 

where, for j = 1, 2, 3 and 4, Bj{ ) denotes the MLE of Bj{-), that is, the unknown parameters are replaced 
by their MLEs. Finally, one can also use bootstrap techniques to obtain directly the biases of fi and (j), 
therefore, denote by Bp_ -^{p.) the vector of 0{n~^) parametric bootstrap biases for fi, and by Bp_ - {(/)) 

the vector of 0{n^^) parametric bootstrap biases for (j) defined in equation (fT3|) . Further, define Bp{fi) 
and Bp{(t)) analogously from the vector Bp{CX) also in equation p3p . Hence, we can also define the 
following bias-corrected estimators of /i and 0: 

Pi = A - Bp- j^ifi) and % = fi - Bp(p), 

and _ _ 

01 = - Bp^^i^) and 02 = - Bpi^^). 



9 



We give in Tables [T] and [2] the most common link functions for gi and g2, respectively, together 
with their first and second derivatives. We believe this will help the practioners that may be interested 
in applying our results. For Table [T] $(•) denotes the standard normal distribution function, f{x) = 
l/-\/27rexp{— l/2x^} is the density of a standard normal distribution and f'{x) — —x/V^exp{—l/2x'^} 
is the derivative of the density of a standard normal distribution. 



Table 1: Values of dfi/drji, cP^/diji for the most common link functions. 



Link function 


Formula 




dji/drji 


Logit 


logW(l - fi)) = 


Vi 




Probit 








Comp. Log-Log 


10g(- 10g(l - ^)) : 


= Vi 


-log(l-^)/(l-M) 


Link function 


Formula 




d^^i/dri'f 


Logit 


logW(l - ^i)) = 


Vi 




Probit 








Comp. Log-Log 


l0g(-lOg(l -/i)) : 


= Vi 


-(l-M)log(l-M)(l + log(l-M)) 



Table 2: Values of d4)/drj2, d'^ip/drj'^ for the most common link functions. 



Link function 


Formula 


d(j)/dr]2 


d'ct>/dr^i 


Identity 




1 





Log 


log((/>) = V2 




<t> 


Square root 


V^=V2 


20 


2 



5 Some special cases 

In this section we consider some special models that commonly arise in the practical use, namely, the 
linear beta regression model, the linear beta regression model with dispersion covariates, the nonlinear 
beta regression model and the nonlinear beta regression model with linear dispersion covariates. Further, 
we study these models in full detail, that is, we give closed form expressions for their score vector, Fisher's 
information matrix and for the ©(n^^) biases of the MLEs of /3, 0, /i and </>. 



5.1 The linear beta regression model 

For linear beta regression models we have, in ([3]), gi{fJ-i) = gifJ-i), where g{-) is some link function, 
<72(0i) = 9^1; and further, we can write Q as 

ff(Mi) ^Vi^ xfP and = 0, 

where (/> > is a constant, i.e., we have that in this case X = X and Z = 1, where X is the matrix of 
covariates with rows given by xf, and the parameters (3 G and (f> G (0,oo). Furthermore, the score 
vector ^ becomes 

Uf3{/3,^) = ^X^Tiy* - ti*), 

where T — Aia,g{d^i / dr^i) and and vt are defined in pagelH Moreover, the matrices P and W 

defined by equations (O and ([7]) becomes, respectively, 



P 



X 
1 
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and 

w — ( 



W=[ ^Ji', Z'' ) , (17) 



with 



= diag I ^^a, ( 



Wf}^ = diag (^{^iiai - - /ii)<^)} 



W^^ = diag(6j). 

Here, Ui and &i are as defined in page [H Therefore, we have the following expression for the Fisher's 
information matrix for the parameter vector C, = (/3"^,</>)'^ 

K{C) = P^WP. (18) 

By means of simple calculations it is possible to conclude that equations ([T6|) and (fT8|) agree with the 
expressions for the score vector and Fisher's information matrix given in Ferrari and Cribari-Neto (2004). 

We now move to the bias correction. Note initially that ^2 given in (|lip vanishes since, for this model, 
both 771 and 772 are linear functions of /? and 9, respectively. Thus, if P and W are the matrices defined 
for this model, equation pip equals simply 

B{c) = {p^wpy^p^w^i, 

which is easy to see that agrees with the finds of Ospina et al. (2006). Further, uj2 in equation ([T^ also 
vanishes which gives us the expression 

c/*(C)-c/(C) + i^V, 

and it is also easy to see that it agrees with the finds of Ospina et al. (2006). Finally, following equation 
P^ . the second order bias of p, can be expressed, in matrix notation, as 

Bifi) = TXB0) + ^SiPfsp, 
with T defined in the begining of this subsection and 6*1 was defined in page [HI 



5.2 The linear beta regression model with dispersion covariates 

This class of models generalizes the linear beta regression models considered in the last subsection by 
letting the precision parameter (j) vary through a linear regression structure. More precisely, for this 
model the equation ([3]) becomes 

51 = Vii = xfP and 52 (0i) = mi = ^I^^, 

where /3 G K*' and 9 eM'^. Then, we have that for this model X — X and Z — Z, where X is the matrix 
of covariates with rows given by xf and Z is the matrix of covariates with rows given by zf . The score 
vector for this model is identical to the one given in expression (O, only that in this case X and Z must 
be replaced by X and Z, respectively. 

Now, the matrix P defined in ^ becomes for this model 

\ Z J ' 

further, let W be the matriz defined in thus, we have the following expression for the Fisher's 
information matrix for the parameter vector ( — {P'^ , 9)^ 

K{C) = P'^WP. 
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Moving to bias correction, since both 771 and 772 arc linear functions of /3 and ^, respectively, we have 
that ^2 in equation (fTT|) vanishes and thus the 0{n^^) bias of the MLEs of the parameter vector ( is 

B{c) = {p^wpy^p^w^i. 

Further, since uj2 in vanishes, the modified score to obtain the estimator through the preventive 
method is given by 

c/*(C) = c/(C)+i^V, 

and, finally, using equations (I14p and psp . the second order biases of fi and 4> are given by 

Bifi) ^ TiXB0) + ^SiPp0, 

and 

5(0) = r2^B(^) + ^S2P9e, 

where Ti and T2 were defined in page HI 5*1 and ^2 were defined in pageO Pf^p and Pee were defined in 
pageini 

5.3 The nonlinear beta regression model 

We now consider the nonlinear beta regression models. For this models, the expressions in ([3]) turns to 
be giifii) = g{fii), where g{-) is some link function, g2{4'i) — 0i, and further, we can write ^ as 

gip-i) = m = f{xf;P) and 0j = 0, 

where > is a constant, i.e., and thus, in this case X remains the same and Z = 1, the parameters 
/3 e M*^ and € (0, 00). Furthermore, the score vector (O becomes 

U^{l3,(f)) = 

where T is as defined in equation p6p . and are defined in page HI Moreover, the matrix P defined by 
equation ([H]), becomes 

' X 
1 

and the matrix W defined by equation ([7]) is actually the same as the matrix given in equation ()17p . 
Therefore, the Fisher's information matrix for the parameter vector C, — {0^ , <j>)'^ can also be written as 
K{C) = P^WP. 

We now turn to bias correction. For this model, the vector ^2 in equation pip can be written as 
^2 = W~^oj2, where 0)2 is given below: 

- ( 'n:f ) ' (19) 

where A^i and N2 are given in equations and ((?7|) of Appendix, and F is given in page [H Thus, 
using P, W and ^2 defined in this Subsection, the 0{n^^) bias for the MLEs of the parameter vector C 
is given by 

B{C) = iP^WP)-'P^W{^,+^2), 

^1 being as defined in formula (jlip . Further, the modified score to obtain the estimator through the 
preventive method is given by 

u*iO = uiO + P^i^i+co2). 

Moreover, the second order bias of /t can be written, following equation (|14p . as 

B{fi) - ^T{2XB0) +F) + ^S^Ppp, 



with 5*1 as defined in pageHl 
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5.4 The nonlinear beta regression model with linear dispersion covariates 

This class of models generalizes the nonlinear beta regression models considered in the last subsection 
by letting the precision parameter (p vary through a linear regression structure, such as the linear beta 
regression with dispersion covariates generalizes the linear beta regression. More precisely, for this model 
the equation ([3]) becomes 

9i{f^t) = Vii = fi^f; f3) and 52(01) = V2t = zfO, 

where /3 S R'' and 6* e M'*. Then, we have that for this model X remaining the same, and Z — Z, where 
Z is the matrix of covariates with rows given by zf . The score vector for this model is identical to the 
one given in expression ([5]), only that, in this case Z must be replaced by Z. 
Now, the matrix P defined in fB]) becomes for this model 

I)' 

further, let W be the matriz defined in ([T]), thus, we have the following expression for the Fisher's 
information matrix for the parameter vector C — {0^ , 9)^ 

K{C) = P^WP. 

Moving to bias correction, for this model the vector ^2 in equation can be written as ^2 = W~^lj2, 
where 0J2 was defined in equation (jl9p . Thus, using P and ^2 defined in this Subsection, the 0{n~^) bias 
for the MLEs of the parameter vector C is given by 

i?(c) = (p^M/p)-ip^M^(a + 6), 

^1 being as defined in formula (fTT|) . Further, the modified score to obtain the estimator through the 
preventive method is given by 

U*{C) = U{C) + P^{u;i+Cj2). 
Moreover, the second order biases of jl and (fi can be written, following equation (|14p . respectively, as 

B{il) = ^Ti{2XB0) +F) + ^SiPpfs, 

and 

B{4>) = T2ZB{e) + ]^S2Pe0. 
with Ti and T2 as defined in pagelH and. Si and S2 as defined in pageO 

6 Numerical results 

In this section we present the results of some Monte Carlo simulation experiments, where we study the 
finite-sample distributions of the MLEs of (3 and 6 along with their corrected versions proposed in this 
paper. The first experiment uses a logit link in a nonlinear model for the regression parameters and a 
log link in a nonlinear model for the precision parameter 

logit^i = /3o + PiXiA + xl^'i, 

log^i = ^0 + ^ia;i,i + xl^^, i = 1, . . . , rt, 

where the true values of the parameters were taken as (3q = 1.5, /3i = 0.5 and (32 = 2; and = 1-7, 6*1 — 
0.7 and 02 — 3. Note also that here the elements of the n x 3 matrix X are: X{f3)i^i — l\X{(3)i^2 — a;i,i, 
and X(/3)i^3 = log(2;2, 1)2^2 V The explanatory variables xi and X2 were generated from the standard 
normal and uniform U(l,2) distributions, respectively, for n — 20,40 and 60. The values of xi and X2 
were held constant throughout the simulations. The total number of Monte Carlo replications was set 
at 5,000 for each sample size. All simulations were performed using the software Ox. 
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Table 3: Simulation results for n — 20. 



Parameter MLE Cox-Snell Firth p-boot np-boot 



f3o 


1, 


.5147 


1.5032 


1, 


.5253 


1.5019 


1, 


.5061 


Bias 


0, 


.0147 


0.0032 


0, 


.0253 


0.0019 


0, 


.0061 


Variance 


0, 


.0628 


0.0627 


0, 


.0667 


0.0627 


0, 


.0677 


MSB 


0, 


.0630 


0.0627 


0, 


.0674 


0.0627 


0, 


.0677 


Pi 


0, 


.5001 


0.5001 


0, 


.5002 


0.5005 


0, 


.5017 


Bias 


0, 


.0001 


0.0001 


0, 


.0002 


0.0005 


0, 


.0017 


Variance 


0, 


.0126 


0.0126 


0, 


.0132 


0.0125 


0, 


.0143 


MSB 


0, 


.0126 


0.0126 


0, 


.0132 


0.0125 


0, 


.0143 


02 


1, 


.9905 


1.9968 


1, 


.9847 


1.9982 


2, 


.0036 


Bias 


-0.0095 


-0.0032 


-0.0153 


-0.0018 


0, 


.0036 


Variance 


0, 


.0119 


0.0118 


0, 


.0129 


0.0117 


0, 


.0125 


MSB 


0, 


.0120 


0.0118 


0, 


.0130 


0.0117 


0, 


.0125 


Oo 


1, 


.9512 


1.7406 


2, 


.1388 


1.7045 


1, 


.8107 


Bias 


0, 


.2512 


0.0406 


0, 


.4388 


0.0045 


0, 


.1107 


Variance 


0, 


.5615 


0.5100 


0, 


.6248 


0.4833 


0, 


.4739 


MSB 


0, 


.6247 


0.5116 


0, 


.8173 


0.4833 


0, 


.4861 


ei 


0, 


.7091 


0.6928 


0, 


.7331 


0.6842 


0, 


.7060 


Bias 


0, 


.0091 


-0.0072 


0, 


.0331 


-0.1584 


0, 


.0060 


Variance 


0, 


.2881 


0.2481 


0, 


.3651 


0.2231 


0, 


.2369 


MSB 


0, 


.2882 


0.2481 


0, 


.3662 


0.2234 


0, 


.2369 


O2 


3 


.0261 


3.0061 


3 


.0440 


3.0043 


3, 


.0757 


Bias 


0, 


.0261 


0.0061 


0, 


.0440 


0.0043 


0, 


.0757 


Variance 


0, 


.0706 


0.0610 


0, 


.0805 


0.0514 


0, 


.0471 


MSB 


0, 


.0720 


0.0611 


0, 


.0825 


0.0514 


0, 


.0528 



In each of the 5, 000 replications, we fitted the model and computed the MLBs /3, 9, its corrected 
versions from the corrective method (Cox and SncU, 1968), preventive method (Firth, 1993) and the 
bootstrap method both of its parametric and nonparametric versions (Bfron, 1979). The number of 
bootstrap replications was set to 600 for both bootstrap methods. 

In order to analyze the results we computed, for each sample size and for each estimator, the mean 
of estimates, bias, variance and mean square error (MSB). Tables [51 1151 present simulation results for 
sample sizes n = 20, 40 and 60, respectively. Moreover, we also considered the estimated values of 
and (f)i, i = 1, . . . ,n, n = 20,40 and 60, along with its corrected estimators presented in Section H) We 
also want to emphasize that, for the bootstrap schemes, two methods to removing the first-order bias 
can be considered: the one induced by the Taylor's expansion and the one given by the bootstrap itself, 
denoted by Si and cti, respectively, where is a surrogate for /i^ and 1/)^. 

Table [3] presents simulation results for sample size n = 20 with respect to the parameters /3 and 6. 
We begin by looking the estimated biases, in absolute value, of the estimators. Initially we note that for 
all parameters the biases of the corrective estimators were smaller than of the original MLBs. However, 
for all parameters the biases of preventive estimators were bigger than of the original MLBs, moreover, 
not only the biases were bigger but also the MSB, which shows that the preventive method does not work 
well for this model, the same phenomenon occurred in Ospina et al. (2006), which collaborates to the 
idea that this method has some problems in beta regression models. We now observe that the MSB of the 
corrective estimators were smaller than those of the MLBs for all parameters, showing that the correction 
is effective. Moving to the bootstrap corrected-estimators we note that the parametric bootstrap had the 
smallest MSB for all parameters, even though the biases were not the smallest. However, the MSB were 
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Table 4: Estimated values of for n = 20. 



i MLE Cox-Snell p-boot p-boot np-boot np-boot 



fij MSExlO^ fij MSExIQS J[^^ MSExlO^ Jt^j MSExlO^ Ji^j MSExlO^ /l^.i MSExlO^ 



1 
i 


u 


.yyo / D 


u. 


.UlUo 


n 

u. 


,yyo ( D 


U, 


,UiUo 


u 


.yyo / / 


u 


ni 07 

.UlU 1 


n 

u. 


,yyo ( ( 


A 
U, 


ni 07 

,U1U ( 


n 
u 


.yyooi 


u 


ni 9s 


A 
U, 


,yyooo 


A 
U, 


ni 99 


o 
It 


n 
u 


.yyoOO 


n 


997fi 
.ZZ 1 


n 

u. 


.yyooi 


A 

U, 


9991 
.ZZZl 


n 
u 


.yyoD4 


u 


.ZZUo 


A 

u. 


.yyooi 


A 
U, 


999"? 


u 




u 


9fi1 fi 
.ZDID 


A 
U, 


.yyooy 


A 
U, 


91^71 


q 
o 


u 


.youoi 


u. 


nnni 


u, 


.youoy 


1 

1, 




u. 


.yOUc5i 




r\AAr\ 

.U44U 


u. 


.youoo 


1 

1, 


r\AQO 

.U4yz 


u. 


.yozoz 




.izoy 


A 

U, 


.youyz 


1 

1, 


1 A'iO 


4 


0, 


.98047 


0, 


,7630 


0, 


.98049 


0, 


.7620 


0, 


.98051 


0, 


.7594 


0, 


98050 


0, 


.7616 


0, 


.98084 


0, 


.8367 


0, 


.98070 


0, 


.8373 


5 


0, 


.95597 


7. 


,5592 


0, 


.95612 


7, 


.5248 





.95621 


7, 


.4940 


0, 


.95609 


7, 


.5313 


0, 


.95755 


8, 


.1211 


0, 


.95635 


8, 


.2285 


6 





.99447 


0, 


,0152 


0, 


.99447 


0, 


.0152 





.99447 





.0152 


0, 


.99448 


0, 


.0151 





.99451 





.0168 


0, 


.99458 


0, 


.0169 


7 


0, 


.97837 


1, 


,1784 


0, 


.97838 


1, 


.1761 


0, 


.97842 


1 


.1737 


0, 


.97838 


1, 


.1783 


0, 


.97897 


1 


.2768 


0, 


.97858 


1, 


.2786 


8 


0, 


.99383 


0, 


,0126 


0, 


.99383 


0, 


.0126 


0, 


.99383 





.0126 


0, 


.99384 


0, 


.0126 





.99387 





.0143 


0, 


.99394 


0, 


.0146 


9 





.99537 


0. 


,0117 


0, 


.99537 


0, 


.0117 





.99537 





.0117 


0, 


.99538 


0, 


.0116 





.99542 





.0133 


0, 


.99547 


0, 


.0129 


10 


0, 


.98111 


1, 


,1470 


0, 


.98116 


1, 


.1389 


0, 


.98120 


1, 


.1357 


0, 


.98116 


1, 


.1418 


0, 


.98183 


1 


.2585 


0, 


.98133 


1, 


.2525 


11 


0, 


.99451 


0, 


,0115 


0, 


.99450 


0, 


.0115 


0, 


.99450 


0, 


.0115 


0, 


.99451 


0, 


.0114 


0, 


.99455 





.0134 


0, 


.99460 


0, 


.0134 


12 





.99002 


0. 


,1023 


0, 


.99002 


0, 


.1022 





.99003 





.1020 


0, 


.99004 


0, 


.1020 





.99013 





.1126 


0, 


.99018 


0, 


.1133 


13 


0, 


.97723 


2, 


,1492 


0, 


.97732 


2, 


.1264 


0, 


.97739 


2, 


.1185 


0, 


.97732 


2, 


.1327 


0, 


.97832 


2, 


.3582 


0, 


.97750 


2, 


.3507 


14 


0, 


.99222 


0. 


,1662 


0, 


.99224 


0, 


.1646 


0, 


.99226 


0, 


.1640 


0, 


.99225 


0, 


.1648 





.99251 


0, 


.1888 


0, 


.99235 


0, 


.1856 


15 





.99169 


0. 


,0522 


0, 


.99167 


0, 


.0523 





.99168 





.0523 


0, 


.99168 


0, 


.0523 





.99178 





.0590 


0, 


.99180 


0, 


.0586 


16 


0, 


.96273 


4, 


,7219 


0, 


.96282 


4, 


.7072 


0, 


.96289 


4, 


.6905 


0, 


.96280 


4, 


.7116 


0, 


.96393 


5, 


.0719 


0, 


.96305 


5, 


.1305 


17 


0, 


.97674 


2, 


,0972 


0, 


.97687 


2, 


.0702 





.97693 


2, 


.0568 


0, 


.97688 


2, 


.0675 





.97758 


2 


.2964 


0, 


.97711 


2, 


.3020 


18 





.95009 


1. 


,3879 


0, 


.95039 


1, 


.3752 





.95056 


1 


.3703 


0, 


.95035 


1, 


.3794 





.95288 


1 


.5040 


0, 


.95061 


1, 


.5142 


19 


0, 


.94657 


1, 


,5431 


0, 


.94695 


1, 


.5243 


0, 


.94712 


1 


.5146 


0, 


.94691 


1, 


.5251 


0, 


.94909 


1 


.6708 


0, 


.94717 


1, 


.6917 


20 


0, 


.92738 


3. 


,6923 


0, 


.92811 


3, 


.6351 





.92840 


3 


.6103 


0, 


.92804 


3, 
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0, 


.93157 


4, 
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0, 


.92828 


4, 


.0608 



Table 5: Estimated values of (p ioi n = 20. 



i 


True 




MLE 


Cox-Snell 


_ P- 


-boot 


np 


-boot 








MSExlO"^ 




MSExlO"^ 


6, ■ 


MSExlO"^ 




MSExlO"^ 


1 


12218 


32123 


48311 


1 9500 




8542.3 


831.41 


13848 


21445 


2 


1512.2 


4977.5 


6716.2 


1567.4 


212.00 


873.67 


97.962 


1991.9 


809 20 


3 


33.710 


55.534 


0308 


35.446 


0082 


33.122 


0.0071 


39.450 


0.0138 


4 


381.74 


652 63 


3.5711 


404.02 


0.8616 


374.10 


6805 


523 56 


2.3280 


5 


49.260 


81.864 


0.0911 


52.560 


0.0236 


49.382 


0.0187 


60.855 


0.0544 


6 


8602.5 


22211 


19727 


8858.7 


1522.2 


6182.9 


308.90 


9874.8 


9677.1 


7 


168.67 


267.63 


0.4576 


174.93 


0.1307 


163.20 


0.1128 


218.12 


0.2947 


8 


5343.4 


11400 


2349.0 


5613.7 


34.533 


4558.3 


121.79 


7036.3 


843.28 


9 


12920 


37410 


90373 


12865 


4081.7 


7660.0 


2134.6 


12609 


58969 


10 


191.66 


343.11 


2.0051 


203.41 


0.4408 


182.27 


0.2683 


248.99 


0.7876 


11 


6293.7 


13636 


3631.7 


6597.4 


0.0518 


5301.0 


180.91 


8269.6 


1213.7 


12 


2235.6 


4604.2 


304.53 


2362.9 


50.021 


2000.5 


23.615 


3035.0 


124.77 


13 


122.77 


237.92 


1.5913 


132.08 


0.2858 


114.87 


0.1369 


152.74 


0.4410 


14 


1288.6 


2744.5 


306.24 


1377.5 


44.117 


1142.9 


15.643 


1843.3 


82.982 


15 


1706.1 


2935.0 


55.657 


1767.5 


13.164 


1589.4 


9.5611 


2382.9 


32.945 


16 


66.492 


106.15 


0.0928 


69.929 


0.0027 


66.068 


0.0234 


83.020 


0.0559 


17 


375.66 


876.02 


33.651 


412.62 


3.7010 


351.16 


3.1617 


520.18 


23.952 


18 


28.969 


57.292 


0.0949 


30.760 


0.0153 


26.626 


0.0068 


30.897 


0.0227 


19 


43.286 


98.112 


1.9018 


50.435 


0.1725 


44.321 


0.1744 


57.730 


2.6463 


20 


24.754 


81.680 


9.3312 


30.351 


0.1859 


22.880 


0.1946 


39.086 


15.773 



Table 6: Simulation results for n — 40. 



Parameter 
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P- 
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1 


.5079 


1, 
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1, 
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0, 


.0337 


0, 
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0, 
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0, 
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0, 
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1, 
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1.9993 


1, 
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2, 


.0014 


2, 


.0069 


Bias 


-0 
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-0.0007 


-0 


.0070 


0, 


.0014 


0, 


.0069 
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.0146 


0.0145 


0, 


.0151 


0, 


.0144 


0, 


.0150 


MSE 


0, 


.0146 


0.0145 


0, 


.0151 


0, 


.0144 


0, 


.0150 


do 


1, 


.7460 


1.6884 


1, 
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1, 


.6972 


1, 
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0, 
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-0.0116 


0, 
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-0 
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0, 
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0, 
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0.1755 


0, 
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0, 
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0, 
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0, 
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0, 
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0, 
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0, 
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Oi 


0, 
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0, 
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0, 
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0, 
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Bias 


0, 


.0264 


0.0110 


0, 


.0480 


0, 
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0, 
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0, 
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0.0696 


0, 


.0816 


0, 


.0697 


0, 


.0727 


MSE 


0, 


.0751 


0.0698 


0, 


.0839 


0, 
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0, 
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02 


3 


.0701 


3.0220 


3 


.1226 


3, 


.0067 


3, 


.0632 


Bias 


0, 


.0701 


0.0220 


0, 


.1226 


0, 


.0067 


0, 


.0632 


Variance 


0, 


.0709 


0.0652 


0, 


.0791 


0, 


.0597 


0, 


.0633 


MSE 


0, 


.0758 


0.0657 


0, 


.0941 


0, 


.0597 


0, 


.0673 



very close to the MSE of the corrective method, and the computation of the parametric bootstrap biases 
is computer intensive, whereas the corrective method is not. Lastly, we observe that for all parameters 
9 the MSE of the nonparametric bootstrap corrected estimators were smaller than those of the MLEs. 
Moreover, for the parameters /3, the MSE of the nonparametric bootstrap corrected estimators were very 
close to those of the MLEs, showing that this method is satisfactory, and is very easy to implement by 
practioners since no parametric assumptions are made. Therefore, for the small sample size n = 20, we 
were able to conclude that corrective method by Cox and Snell (1969) was successfully applied, as well 
as the bootstrap correction. 

TableSlpresents the simulation results for sample size n = 20 with respect to the parameter /i. We did 
not consider the preventive estimators since its estimation was poor. The MSE of both nonparametric 
bootstrap correction schemes are the biggest ones for all parameters. The best estimators, with respect 
to MSE, are the parametric bootstrap schemes. Finally, the MSE of the corrective method are always 
least or equal those of the MLE. Note also that, for the for the means, all the values are not distant 
from each other. These results show that, even for small sample, we are able to obtain a satisfactory 
bias correction for both the corrective method and the parametric bootstrap method. 

Table [5] presents the simulation results for sample size n — 20 with respect to the parameter 4>. The 
bootstrap correction scheme based on Taylor's expasion is not presented in this Table, since its perfor- 
mance was worse than those of the regular bootstrap correction, for both parametric and nonparametric 
schemes. The first remarkable fact is that if the true value of is large, the bias of the corrective 
method is small when compared with the others estimators, and the bias of the parametric bootstrap 
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Table 7: Estimated values of /x for n = 40 for observations 1 until 20. 
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Table 8: Estimated values of /x for n = 40 for observations 21 until 40. 
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Table 9: Estimated values of <^ for n = 40 for observations 1 until 20. 
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Table 10: Estimated values of <^ for n = 40 for observations 21 until 40. 
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Table 11: Simulation results for n — 60. 



Parameter MLE Cox-Snell Firth p-boot np-boot 



f3o 
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Variance 
MSE 


1.5046 
0.0046 
0.0156 
0.0156 


1.5025 
0.0025 
0.0156 
0.0156 


1.5063 
0.0063 
0.0159 
0.0159 


1.5016 
0.0016 
0.0156 
0.0156 


1.5005 
0.0005 
0.0160 
0.0160 


/3i 
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MSE 


0.4995 
-0.0005 
0.0030 
0.0030 


0.4996 
-0.0004 
0.0030 
0.0030 
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-0.0004 
0.0031 
0.0031 
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-0.0003 
0.0030 
0.0030 
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0.0000 
0.0033 
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0.0037 
0.0037 
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-0.0061 
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0.0956 
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0.0066 
0.0993 
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0.6966 
-0.0034 
0.0589 
0.0589 


0.6974 
-0.0026 
0.0551 
0.0551 


0.6955 
-0.0045 
0.0650 
0.0650 


0.6988 
-0.0012 
0.0543 
0.0543 
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0.0015 
0.0581 
0.0581 


O2 
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Variance 
MSE 


3.0630 
0.0630 
0.0265 
0.0305 


3.0107 
0.0107 
0.0260 
0.0262 


3.1184 
0.1184 
0.0274 
0.0414 


2.9955 
-0.0045 
0.0257 
0.0257 


3.0073 
0.0073 
0.0271 
0.0271 



scheme is considerably larger than the ones of the corrective method. If the true value of cj) is small, the 
parametric bootstrap outperforms the corrective method with respect to bias. The maximum likelihood 
estimator had the worst perfomance with respect to bias and MSE. The parametric bootstrap had, in 
general, the best performance with respect to MSE. Therefore, we conclude that the corrective method 
was satisfactory, but overall, the parametric bootstrap had the best performance. It is also noteworthy 
that, for the precision parameter 0, the correction schemes worked very well, and, therefore, their use 
produces an improved estimation. 

In Table El the results for sample size n = 40 are presented with respect to the parameters (3 and 
6. As observed in Table [3l the corrective method and the parametric bootstrap scheme had the best 
performance with respect to both bias and MSE. The nonparametric bootstrap scheme had performance 
worse than the MLE for the parameters /3, and better than the MLE for the parameters 9, with respect 
to both bias and MSE. For this sample size, we observe again that the preventive estimator had a poor 
performance. The results show that both the corrective method and the parametric bootstrap produce 
better estimators than the MLE, with respect to bias and MSE, and again, that the estimates produced 
by the preventive method are not satisfactory. Finally, the behavior of the nonparametric bootstrap 
estimator is the same as the one for sample size n = 20. 

Tables [7] and [S] contains the results for sample size n = 40 with respect to the parameter ^. Initially 
we note that, based on MSE, the corrective estimators had the best performance. We also observe that 
the maximum likelihood and parametric bootstrap estimators had similar performance. Finally, the 
nonparametric estimators had the worst performance with regard to MSE. 

Tables 1^ and [TUl sumarize the results for sample size n = 40 with respect to the parameter (p. For this 
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Table 12: Estimated values of for n = 60 for observations 1 until 30. 
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Table 13: Estimated values of fi iov n = 60 for observations 31 until 60. 



i MLE Cox-Snell p-boot p-boot np-boot np-boot 



A/[SExl05 fl, MSExIQS 71^^ A/ISExIQS p^^^ MSExlO^ /i^ , MSExlQ-'^ Jl^ ,^ A/TSExlO^ 



01 


n 
u 


QQ/1 1 K 


u 


.Ulo4 


A 
U. 


QQ/1 1 


A 
U. 


ni s/i 



u 


QQ/1 1 ^ 

.yy4io 


A 
U 


ni 


A 
U. 


QQ4 1 

.yy4io 


A 
U, 


ni Rf; 


A 
U 


QQ/1 1 7 

,yy4i i 


u 


ni Q7 
.uiy i 


A 
U. 


QQ/1 1 R 

.yy4io 


A 
U. 


ni Q7 
,uiy / 




u 


QQ91 R 


u 


.uooy 


n 


QQ91 7 

.yyzi / 


A 

U. 


.uoo / 



u 


QQ91 7 

.yyzi ( 




u 


.U00c5 


A 

U. 


0091 7 

.yyzi / 


A 

U, 


.UOOo 




u 


QQOOA 

,yyzzu 



u 


.uoyi 




u, 


0091 Q 

.yyzio 




u. 


,uoyo 


00 


u 


07700 
.y M 00 


u 




A 

u. 


Q77qQ 

.y ( ( oy 




u. 


.Z ( 00 


n 
u 


Q77qQ 
.y M 00 


A 

u 


97/1 1 
.Z / 4i 




u. 


077qa 
.y ( ( 00 


A 
U, 


97/19 
,Z ( 4Z 


A 

u 


Q77/I A 

,y M 4u 



u 


.ZoOo 


A 

u. 


077q0 

,y / / oy 


A 

u. 


98/19 
,Zo4Z 




n 
u 


.yooDo 


u 


. IZDo 


A 

u. 


.yoooy 


A 

u. 


1 97n 

.iZ / u 


A 

u 


.yooDo 


n 
u 


1 979 
.IZ 1 Z 


A 

u. 


.yoooo 


A 
U, 


1 979 
, IZ J Z 


A 

u 


QOqCQ 

,yoooy 




u 


1 9SQ 

. izoy 


A 

u. 


QcqcQ 

,yoooy 


A 

u. 


1 9Qn 
,izyu 


00 


u. 


.y ( uio 


1 
± 


.Z ( 00 


n 

u. 


.y /uzD 


1 

1, 


.ZOOO 



u 


.y (UZD 


1 
± 


97nt^ 
.Z ( uo 


A 

U, 


07094 

.y iuz4 


1 

1, 


971 K 
.Z / 10 




u 


Q7Aqo 

.y ( uoz 


1 
± 


.ouoo 




u, 


07091 

,y /uzo 


1 


11 96 


OD 


u 


Q/l f^Q'? 

.y4oyo 


q 



91 nK 


A 

u. 


.y4DUD 


q 


9nf;Q 


n 
u 


Q/l KQ/I 
.y4DU4 


q 



91 '?Q 

.zioy 




u. 


0/1609 

.y4ouz 


q 


91 

,ZiOO 


A 

u 


Q/l KQQ 

.y4ouy 


q 



11 Q9 

.oiyz 


A 

u. 


Q/l p;Q8 

,y4oyo 


q 

0. 


119Q 

,oozy 


i 


u 


.yoyoi 




'?fi9/l 


A 

u. 


.yoyoo 



z. 


.0010 


u 


.yoyoo 



z 


.000 t 


A 

u. 


OtcQRq 
.yoyoo 


9 
z. 


,00 / 4 


A 

u 


QCQ79 

.yoy i z 


9 
z 


A 1 n7 

.41U t 


A 

u. 


QCQGO 

,yoyou 


9 

z. 


,4ZZ0 


00 


u. 


Q009Q 


a 

0. 


.DDUy 


n 


0Q9f;i 
.yozoi 


Q 


Kn77 




u. 


.yozDi 


Q 

0. 


.DoUi 


A 

u, 


0q9c:7 

.yozo / 


(5, 


RIRR 




u. 


01901 

.yozyi 




y. 


m94 

.0UZ4 




u, 


,yozoi 




y, 


1791 
,0 / ZO 




u 


QO/I 07 





■^9/1 


A 

u. 


OQ/I/IR 
.yo44D 


r 
0. 


1 7Q 


u 


Q'?/l/l/l 
.yo444 





.ozyo 




u. 


Q'^/l/l 1 
.y044i 


0, 


,00iy 




u 


Oq/I/IQ 

.yo44y 





A 79S 
.4 / Zo 


A 

u. 


Oq/I qc 
,y0400 


0. 


/IQ/18 

,4y4o 


/in 


u 


.youo i 


1 
± 


"i^i n 

.oOlU 


A 

u. 


QCAQ/I 

.youy4 


1 

1 , 


.ooUo 


n 
u 


Qf;nQ9 
.youyz 


1 
i 


.oOoZ 


A 

u. 


0RnQ9 

.youyz 


1 

1 , 


,0000 


A 

u 


Q6AQC 

.youyo 


1 
± 


1Q1 9 

.oyiz 


A 

u. 


QGAQA 

,youyu 


1 

1 , 


1Qi^6 

,oyoo 


41 


u 


.yc5iui 


U. 


1 Q71 

.iy ( i 


n 


001 09 

.yoiuz 


A 

U, 


1 Q79 

.ly / Z 




u. 


081 n9 
.y<5iuz 




u. 


1 07'=; 
.iy ( 


A 

u, 


OQI 09 

.yoiuz 


A 

U, 


1 Q7K 

,iy /D 




u. 


081 m 
.y<5iuo 




u. 


9nni 
.zuuo 




u, 


oai 01 
,yoiuo 




u, 


9nn'^ 
,zuuo 


/1 9 


u 


QS91 

.yozio 


u 


1 '?sf; 

. looD 


A 

u. 


Q891 

.yozio 




u. 




u 


QS91 


A 

u 


1 '?Qn 




u. 


Q891 
.yoZiO 


A 
U, 


1 '^Ql 

,ioyi 




u 


QS91 A 

.yozi4 




u 


1 41 n 

. 141U 




u. 


QS91 A 

,yozi4 


A 

u. 


1/111 


/I 


u 


QS/I "IQ 

.yo4oy 


u 


1 197 


A 

u. 


QS/i/in 

.yo44U 


A 

u. 


1 "^97 
.LoZ i 


A 

u 


QO/I /I A 
.yo44U 


A 

u 


1 '?9Q 

.lozy 


A 

u. 


QO/I qo 

.yo4oy 


A 
U, 


,iOOU 




u 


QS/I/1 1 

.yo44i 




u 


1 ic:i 
.1000 


A 

u. 


QS/I/1 1 

,yo44i 


A 

u. 


1 if^f; 
,1000 


44 
44 


u. 


QQCQI 

.y<50<5i 


u. 


.UDo4 


n 


ooc;qo 
.yoOou 


A 

u, 


.UDoO 




u. 


QQC7Q 

.yc50 ( y 




u. 


.UDoD 


A 

u, 


oaK7o 
.yoo / y 


A 

U, 


,UD00 




u. 


QQCQI 

.y<50c5i 




u. 


.UDOZ 




u, 


,yoOoi 




u, 


,UDOO 


/I ^ 


u 


.youOo 


1 

1 


.UD < z 


A 

u. 


.yoD ( u 


1 

i. 


.UDOO 


u 


.yoooy 


1 

i 


.uo / 


A 

u 


.yoooo 


1 

i. 


n6Qi 

,UOoi 


n 
u 


Q6671 

.yoo ( i 


i 


nQ6/i 

.Uo04 


A 

u 


,yooo ( 


1 

i. 


n80f^ 
,uoyo 


/If; 

40 


n 
u 


.y 1 ozo 


n 
u 


/I ^9fi 
.40ZD 


A 

u. 


Q7i^'J1 

.y / ooi 


A 
U. 


A f;9/i 

.40Z4 


A 

u 


Q7C qn 

.y t oou 


A 

u 


.40oZ 


A 

u. 


Q7c;oo 

.y / ozy 


A 
U, 


,4004 


A 

u 


Q7C q9 

.y i ooz 


A 

u 


/I fin/I 

.40U4 


A 

u. 


07cqA 

,y / oou 


A 

u. 


4611 
,4010 


47 
4 ( 


n 
U 


QQnn4 
.yyuu4 


u. 


n994 
.UZZ4 


n 


.yyuuo 


A 

U, 


0994 
.UZZ4 




u. 


QQAAO 

.yyuuo 




u. 


n994 
.UZZ4 


A 

u, 


.yyuuo 


A 

U, 


0994 
,UZZ4 




u. 


.yyuu4 




u. 


n99a 

.UZZc5 




u. 


00004 

,yyuu4 




u, 


0990 

,uzzy 


/IS 
4o 


n 
u 


QO/l 7/1 

.yo4 / 4 


u 


nsn9 

.UoUZ 


A 

u. 


OQ/1 7/1 
.yo4 ( 4 




u. 




u 


QO/I 7q 
.yo4 / 


A 

u 


nsn/i 

.UoU4 




u. 


QQ/1 7q 
.yo4 ( 


A 
U, 


nsn4 

,UoU4 




u 


QO/I 7/1 

.yo4 / 4 




u 


nsi 7 




u. 


084 7K 

,yo4 ( 


A 

u. 


n8i 7 

,Uol ( 


49 





.99354 


0, 


.0221 


0, 


.99354 


0, 


.0221 





.99354 


0, 


.0221 


0, 


.99354 


0, 


,0221 





99356 

1 t/ t/ t_f KJ \J 





.0235 


0, 


99355 


0, 


,0235 


50 





.98988 





.1142 


0, 


.98989 


0, 


.1138 





.98990 





.1139 


0, 


.98990 


0, 


,1140 





.98995 





.1220 


0, 


,98992 


0, 


,1217 


51 


0, 


.98646 


0, 


.0657 


0, 


.98645 


0, 


.0658 


0, 


.98645 


0, 


.0659 


0, 


.98645 


0, 


,0659 


0, 


.98646 


0, 


.0667 


0, 


,98647 


0, 


,0667 


52 


0, 


.92288 


7, 


.7763 


0, 


.92312 


7, 


.7688 


0, 


.92309 


7, 


.7860 


0, 


.92306 


7, 


,7895 





.92316 


8, 


.0112 


0, 


,92297 


8, 


,0444 


53 





.94860 


3 


.0225 


0, 


.94873 


3 


.0171 





.94872 


3 


.0242 


0, 


.94870 


3 


,0258 





.94880 


3 


.1809 


0, 


,94867 


3 


,1969 


54 


0, 


.98109 


0, 


.2037 


0, 


.98110 


0, 


.2038 


0, 


.98109 


0, 


,2042 


0, 


.98109 


0, 


,2043 


0, 


,98111 


0, 


.2071 


0, 


,98110 


0, 


,2074 


55 





.95698 


1 


.8418 


0, 


.95708 


1, 


.8398 





.95706 


1 


.8435 


0, 


.95705 


1, 


,8444 





,95709 


1 


.8839 


0, 


,95702 


1, 


,8899 


56 





.94157 


4, 


.1167 


0, 


.94173 


4, 


.1093 





.94171 


4, 


.1190 


0, 


.94169 


4, 


,1212 





.94180 


4, 


.3166 


0, 


,94164 


4, 


,3388 


57 


0, 


.95495 


2, 


.0132 


0, 


.95505 


2, 


.0114 


0, 


.95503 


2, 


.0156 


0, 


,95502 


2, 


,0165 





.95506 


2, 


.0675 


0, 


,95499 


2, 


,0745 


58 


0, 


.99070 


0, 


.0199 


0, 


.99069 


0, 


.0199 


0, 


.99069 


0, 


.0200 


0, 


,99069 


0, 


,0200 





.99070 





.0205 


0, 


,99071 


0, 


,0206 


59 


0, 


.98021 


0, 


.2164 


0, 


.98021 


0, 


.2166 


0, 


.98021 





.2171 


0, 


.98021 


0, 


,2172 





.98024 


0, 


.2291 


0, 


,98022 


0, 


,2294 


60 





.96732 





.8075 


0, 


.96736 


0, 


.8076 





.96735 





.8092 


0, 


.96734 


0, 


,8096 





.96737 





.8327 


0, 


,96734 


0, 


,8349 



Table 14: Estimated values of (f) for n = 60 for observations 1 until 30. 
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Table 15: Estimated values of <j)ior n = 60 for observations 31 until 60. 
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sample size the MLE had the biggest bias and MSE, and hence was the poorest estimator considered in 
this table. On the opposite direction, the parametric bootstrap estimator had the smallest MSE, followed 
by the corrective method. However, if the true value is large, the corrective method has bias smaller 
than the parametric bootstrap. It is also noteworthy that for large values of </> the parametric bootstrap 
method tends to underestimate the parameter. We note, again, that both the corrective and parametric 
bootstrap estimators provide better results than the MLE, with respect to both bias and MSE. 

Table [TT] summarizes the results for sample size n = 60 with respect to the parameters /3 and 9. We 
begin by noting that even for a large sample size the preventive estimator had the poorest performance. 
The parametric bootstrap estimator had, in general, a better performance than the MLE, with respect 
to both bias and MSE. The nonparametric bootstrap estimator showed the same behavior it did for the 
sample sizes n — 2Q and n = 40, i.e., had better perfomance for 9 and worse performance for /3 when 
compared to the MLE. Finally, the best performance was obtained by the corrective estimator, which 
had, in general, the smallest bias and MSE. Hence, the corrective method can be effective even for large 
sample sizes. 

Tables [12] and [13] present the results for sample size n = 60 with respect to the parameter /i. We 
note that both nonparametric estimators had the worst performance with respect to MSE. Further, the 
corrective method presented the best performance with respect to MSE, nevertheless it was not much 
better than the MLE. The parametric bootstrap had a performance worse than both maximum likelihood 
and the corrective estimator, but better than the nonparametric estimator. But, overall, all estimators 
were similar. 

Tables [T4l and [T5l present the results for sample size n = 60 with respect to the parameter 0. The MLE 
had the worst performance with respect to both bias and MSE. Considering the MSE, the best estimator 
was the parametric bootstrap, which had consistently the smallest MSE. The parametric bootstrap was 
closely followed by the corrective estimator which had a similar performance, except when the true 
value of 4) was large, in this case, the MSE of the parametric bootstrap was considerably smaller than 
the MSE of the corrective method, but, on the other hand, the bias of the corrective method was also 
considerably smaller than the bias of the parametric bootstrap. Finally, the nonparametric bootstrap 
outperformed the MLE. Therefore, we conclude that even for sample size large the corrective method 
and the parametric bootstrap may produce improved estimators. 

We have performed a second set of simulations, but now we have two goals. First, to check the 
performance of the MLE against the proposed estimators, but secondly, we will consider a model which 
is linearizable, and we will compare the behavior of the estimators for both the linear and non-linear 
fits. We hope this will give another motivation for the usage of the class of nonlinear regression models 
that we are introducing here in this paper. This strategy was also used by Cook et al. (1986). In this 
experiment we consider a linearizable nonlinear model with logit link for the regression parameters and 
a linearizable nonlinear model with identity link for the precision parameter: 

logit^j = /?oa;f ' , 

(j)i=9(ixl\ i = l,...,n, 

where the true values of the parameters were taken as /?o — 0.7 and /3i — 0.5; and 9q = 100 and 9i — 2. 
Note also that here the elements of the n x 2 matrix X are: X{l3)i^i = xf^ and X{(3)i^2 = Po '^og{xi){xf^). 
The explanatory variable x was generated from random draws of an exponential distributions with mean 
1 for n = 20, 40 and 60. The values of x were held constant throughout the simulations and were the 
same for both the linearizable and the nonlinear model. The total number of Monte Carlo replications 
was set at 5,000 for each sample size. The number of bootstrap replications was set to 600 for both 
bootstrap methods. 

To fit the linearized model, we consider the model in the following sense: 

gif^i) = 70 + PiZi, 
h{(f>i) = Co + 9iZi, i = l,...,n, 

where 70 = log(/3o), Co = log(6'o), Zi = log(a;i), g{-) = log(logit(-)) and finally, h{-) = log(-). Hence, we 
can obtain the MLE of Po in the linearized model simply by exponentiation of 70, and analogously to 
obtain the MLE of 6*0. 
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Table 16: Simulation results for n — 20. 



Model Parameter MLE Cox-Snell Firth p-boot np-boot 



Linear 


Bias 

Variance 
MSE 


0.7008 
0.0008 
0.0039 
0.0039 


0.7020 
0.0020 
0.0039 
0.0039 


0.6998 
-0.0002 
0.0039 
0.0039 


0.7022 
0.0022 
0.0039 
0.0039 


0.7040 
0.0040 
0.0040 
0.0040 




Pi 
Bias 

Variance 
MSE 


0.5025 
0.0025 
0.0115 
0.0116 


0.5014 
0.0014 
0.0115 
0.0115 


0.5035 
0.0035 
0.0116 
0.0116 


0.5014 
0.0014 
0.0115 
0.0115 


0.4991 
-0.0009 
0.0132 
0.0132 




Oo 

Bias 

Variance 
MSE 


139.37 
39.370 
3849.4 
5399.1 


111.43 
11.430 
2444.3 
2575.0 


171.78 
71.780 
5917.0 
11070 


107.73 
7.7300 
2307.3 
2367.1 


109.10 
9.1000 
2433.6 
2516.3 




ei 

Bias 

Variance 
MSE 


2.0871 
0.0871 
0.1548 
0.1624 


2.0215 
0.0215 
0.1489 
0.1494 


2.1603 
0.1603 
0.1630 
0.1887 


1.9945 
-0.0055 
0.1488 
0.1488 


1.9791 
-0.0209 
0.1530 
0.1535 


Nonlinear 


Po 
Bias 

Variance 
MSE 


0.7016 
0.0016 
0.0038 
0.0038 


0.7006 
0.0006 
0.0038 
0.0038 


0.7028 
0.0028 
0.0038 
0.0038 


0.7005 
0.0005 
0.0038 
0.0038 


0.7014 
0.0014 
0.00394 
0.0039 




Pi 

Bias 

Variance 
MSE 


0.5007 
0.0007 
0.0115 
0.0115 


0.4997 
-0.0003 
0.0115 
0.0115 


0.5008 
0.0008 
0.0115 
0.0115 


0.4998 
-0.0002 
0.0115 
0.0115 


0.4977 
-0.0023 
0.0128 
0.0129 




Oo 

Bias 

Variance 
MSE 


140.33 
40.330 
3822.6 
5448.8 


101.14 
1.1400 
1961.5 
1962.8 


199.46 
99.460 
7906.2 
17798 


83.136 
-16.864 
1310.9 
1595.3 


79.011 
-20.989 
3575.2 
4015.7 




Oi 

Bias 

Variance 
MSE 


2.0923 
0.0923 
0.1558 
0.1643 


2.0264 
0.0264 
0.1496 
0.1503 


2.2071 
0.2071 
0.1667 
0.2096 


1.9991 
-0.0009 
0.1482 
0.1482 


1.9863 
-0.0137 
0.1536 
0.1538 



Table [TBI presents the results for n = 20. The best estimator, with respect to both bias and MSE, 
was the parametric bootstrap. The preventive estimator had the poorest performance. Considering both 
the linear and nonlinear models, the nonparametric bootstrap had better performance than the MLE 
for the parameter 9 and worse for the parameter (3. The corrective estimator had a good performance, 
in the sense that only the parametric bootstrap had outperformed it. Now, comparing the linear and 
nonlinear models, we note that the MLE was similar for both models. Comparing the corrective and 
parametric bootstrap estimators for both models, the estimators from the nonlinear model had the best 
performance. The nonparametric bootstrap estimator had a similar performance in both models, except 
for ^0 which was considerably worse in the nonlinear model. Hence, we conclude that for the estimators 
of interest, i.e., corrective and parametric bootstrap, the nonlinear model had the best performance and 
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Tabic 17: Simulation results for n — 40. 



Model 


Parameter 


MLE 


Cox-Snell 


Firth 


p-boot 


np-boot 


Linear 


Po 
Bias 

Variance 
MSE 


0.7001 
0.0001 
0.0023 
0.0023 


0.7013 
0.0013 
0.0023 
0.0023 


0.6990 
-0.0010 
0.0023 
0.0023 


0.7018 
0.0018 
0.0023 
0.0023 


0.7015 
0.0015 
0.0023 
0.0023 




Pi 
Bias 

Variance 
MSE 


0.5015 
0.0015 
0.0097 
0.0097 


0.4999 
-0.0001 
0.0096 
0.0096 


0.5030 
0.0030 
0.0097 
0.0097 


0.4994 
-0.0006 
0.0097 
0.0097 


0.4996 
-0.0004 
0.0098 
0.0098 




Bias 

Variance 
MSE 


116.87 
16.870 
1081.2 
1365.7 


104.55 
4.5500 
861.77 
882.46 


130.11 
30.110 
1348.1 
2254.8 


104.18 
4.1800 
882.60 
900.09 


104.25 
4.2500 
862.99 
881.06 




Oi 

Bias 

Variance 
MSE 


2.0346 
0.0346 
0.0663 
0.0675 


2.0076 
0.0076 
0.0647 
0.0647 


2.0627 
0.0627 
0.0682 
0.0721 


2.0029 
0.0029 
0.0662 
0.0662 


1.9939 
-0.0061 
0.0657 
0.0657 


Nonlinear 


Po 
Bias 

Variance 
MSE 


0.7002 
0.0002 
0.0023 
0.0023 


0.6999 
-0.0001 
0.0023 
0.0023 


0.7005 
0.0005 
0.0023 
0.0023 


0.6999 
-0.0001 
0.0023 
0.0023 


0.7000 
0.0000 
0.0023 
0.0023 




Pi 

Bias 

Variance 
MSE 


0.5017 
0.0017 
0.0096 
0.0096 


0.5001 
0.0001 
0.0095 
0.0095 


0.5029 
0.0029 
0.0096 
0.0096 


0.5002 
0.0002 
0.0095 
0.0095 


0.4998 
-0.0002 
0.0097 
0.0097 




Oo 

Bias 

Variance 
MSE 


116.85 
16.850 
1057.6 
1341.7 


100.38 
0.3800 
775.37 
775.51 


137.00 
37.000 
1466.5 
2835.2 


97.198 
-2.8020 
724.96 
732.81 


97.435 
-2.5650 
766.40 
772.98 




Oi 

Bias 

Variance 
MSE 


2.0327 
0.0327 
0.0631 
0.0642 


2.0054 
0.0054 
0.0616 
0.0616 


2.0726 
0.0726 
0.0649 
0.0702 


1.9997 
-0.0003 
0.0617 
0.0617 


1.9929 
-0.0071 
0.0627 
0.0627 



should be used instead of its linearized version if one aims to obtain better performance with respect to 
MSE and bias. 

Table [17] presents the results for rt — 40. The best performance, with respect to both bias and 
MSE, was achieved by the corrective estimator. Considering both the linear and nonlinear models, 
the nonparametric bootstrap had better performance than the MLE for the parameter 9 and worse 
for the parameter (3. The preventive estimator had the worst performance. The parametric bootstrap 
estimator had a good performance, in the sense that only the corrective estimator had a best performance. 
Now, comparing the linear and nonlinear models, we note that the MLE was similar for both models. 
Comparing the corrective, parametric bootstrap and nonparametric estimators for both models, the 
estimators from the nonlinear model had the best performance. Here, we conclude that, for this sample 
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Table 18: Simulation results for n — 60. 



Model 



Parameter MLE Cox-Snell Firth p-boot np-boot 



Linear 



/3o 
Bias 

Variance 
MSB 



0.6718 
-0.0282 
0.0101 
0.0109 



0.6768 
-0.3232 
0.0076 
0.0082 



0.6712 
-0.0288 
0.0101 
0.0109 



0.7322 
0.0322 
0.0109 
0.0119 



0.6659 
-0.0341 
0.0146 
0.0158 



/3i 
Bias 

Variance 
MSB 



0.5435 
0.0435 
0.0249 
0.0268 



0.5327 
0.0327 
0.0168 
0.0178 



0.5442 
0.0442 
0.0249 
0.0268 



0.4705 
-0.0295 
0.0200 
0.0209 



0.5582 
0.0582 
0.0422 
0.0455 



Co 

Bias 

Variance 
MSB 



100.99 
0.9900 
1471.2 
1472.2 



93.896 
6.1040 
1270.5 
1307.7 



108.46 
8.4600 
1703.6 
1775.2 



108.93 
8.9300 
1307.1 
1386.8 



93.676 
-6.3240 
1299.8 
1339.8 



Bias 

Variance 
MSB 



1.9216 
-0.0784 
0.1029 
0.1090 



1.9064 
-0.0936 
0.1032 
0.1120 



1.9358 
-0.0642 
0.1059 
0.1100 



2.0261 
0.0261 
0.1629 
0.1636 



1.8697 
-0.1303 
0.1693 
0.1862 



Nonlinear /So 0.6977 0.6975 0.6981 0.6976 0.6973 

Bias -0.0023 -0.0025 -0.0019 -0.0024 -0.0027 

Variance 0.0033 0.0033 0.0033 0.0033 0.0035 

MSB 0.0033 0.0033 0.0033 0.0033 0.0035 



Pi 
Bias 

Variance 
MSB 



0.4818 
-0.0182 
0.0212 
0.0212 



0.4925 
-0.0075 
0.0058 
0.0059 



0.4716 
-0.0284 
0.0415 
0.0416 



0.4826 
-0.0174 
0.0201 
0.0201 



0.4862 
-0.0138 
0.01867 
0.01869 



Co 

Bias 

Variance 
MSB 



109.19 
9.1900 
635.41 
719.92 



100.93 
0.9300 
451.46 
452.33 



120.89 
20.890 
777.70 
1214.2 



99.256 
-0.7440 
532.02 
532.57 



97.975 
-2.0250 
522.11 
526.21 



Bias 

Variance 
MSB 



l.c 
-0.0011 
0.0323 
0.0323 



1.9839 
-0.0161 
0.0320 
0.0323 



2.0256 
0.0256 
0.0373 
0.0380 



1.9936 
-0.0064 
0.0333 
0.0333 



1.9776 
-0.0224 
0.0317 
0.0322 



size, the corrective estimator along with the bootstrap based estimators had a better performance when 
modelling the nonlinear model, than when modelling its linearized version. 

Table [18] presents the results for n = 60. The best performance, with respect to both bias and 
MSB, was achieved by the corrective and parametric bootstrap estimators, however, their estimates 
were not far from the ones obtained by the MLB, this is probably due to the large sample size. The 
preventive estimator had the worst performance and the nonparametric bootstrap estimator had a similar 
performance to the MLB. It is noteworthy that the best estimates were obtained by the corrective 
estimator in the nonlinear model, followed by the parametric bootstrap estimator also in the nonlinear 
model. Comparing the linear and nonlinear models, we note that the MLB was similar for both models. 
Comparing the corrective, parametric bootstrap and nonparametric estimators for both models, the 
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estimators from the nonlinear model had the best performance. Therefore, we conclude that for all 
sample sizes the nonlinear model should be preferred to its linearized version. 

7 Application to real data 

We now present an application of the linear beta regression model with dispersion covariates, the par- 
ticular case considered in Section [3 The source of the data is Prater (1956). We want to model the 
proportion of crude oil that is converted to petrol after fractional distillation. We have two explanatory 
variables for this model. The first is the level of crude oil, where 10 different possible levels correspond 
to the proportion of crude oil that was vaporized. The second is the temperature in Fahrenheit at which 
all petrol that is contained in the amount of crude oil vaporizes. Ferrari and Cribari-Neto (2004) and 
Ospina et al. (2006) used this data as an illustration of the linear beta regression model and of some 
bias-correction schemes for the model, respectively. 

The sample size is n = 32. The model specification consists of two parts as seen in equation ([3]). 
The first, which is related to the mean, includes an intercept (xi — 1), 9 different dummy variables 
{x2, ■ ■ ■ ,xio) to represent the 10 possible different situations for the level of crude oil and the covariate 
xii, measuring the temperature in Fahrenheit degrees at which all petrol vaporizes. The second, which 
is related to the precision parameter, includes an intercept (zi = 1) and the covariate Z2 ~ xn. 

Table 19: Statistics and p-values of the LRT and ST. 



Test 


Stat, value 


p- value 


LRT 


4.35902 


0.03681 


ST 


6.57124 


0.01036 



The logit link function was used to relate the mean of the response variable to the linear predictor, 
and the log link function was used to relate the precision parameter to its linear predictor. The unknown 
coefficients were estimated through maximum likelihood using the quasi-Newton optimization method 
BFGS (see, for instance. Press et al., 1992) with analytical derivatives. The corrective (based on Cox 
and Snell, 1968), preventive (based on Firth, 1993) and bootstrap bias corrected (based on Efron, 1979) 
bias corrected schemes considered in Sections [3] and H] were also computed. 

We can then write the model we are considering as 

logit(/ii) = /3i + (32X21 + /^sccsj + PiXii + P^x5i + Pexet + Prxji + Psxsi + Pgxgi + /3ioa;ioi + Piixm, 

log{(j)i) ^ 9i + 92Z21, i = l,...,32, 

where logit(a;) = log(a;/(l — x)). 

Looking at Table [50] we see that 6*2 has an standard error of 0.00361. Thus, if we use the asymptotic 
normality to test if 62 = 0, one would obtain the p- value 5.65 x 10"^, which indicates the significance 
of this parameter. As a further study on this direction we considered the likelihood ratio test (LRT) 
and the score test (ST), to test the null hypothesis that the true model does not contain 62 (that is, 
02 = 0), which was the model considered by Ferrari and Cribari-Neto (2004), and Ospina et al. (2006), 
versus, the full model considered above (that is, 6*2 7^ 0). The value of the statistics together with the 
p- value obtained by the comparison with the quantiles of a Xi distribution are given below in Table 1191 
By looking at Table [H] we may conclude that the null hypothesis (that 6*2 — 0) should be rejected, thus 
showing that the general model introduced in Section 2 may be useful to practioners. 

Table [201 presents the maximum likelihood estimates along with their corrected versions and the 
corresponding standard errors. The adjusted versions are the corrective bias-corrected estimator (Cox- 
Snell), the preventive bias-corrected estimator (Firth), the parametric bootstrap bias-corrective estimator 
(p-boot) and the nonparametric bootstrap bias-corrected estimator (np-boot). It can be seen that the 
estimates and corrected estimates for the parameters /3's are very similar, however, for the parameters 
0's, we observe some difference between the maximum likelihood estimate and the corrected ones. 

In Table [51] we give the maximum likelihood estimates of fi and (j) together with their corrected 
versions. We then see that there are not large differences between the estimates of fi and the corrected 
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Table 20: Estimated vahics of the parameters and their standard errors. 



Parameter 


MLE 


Cox-Snell 


Firth 


p-boot 


np-boot 


Ml 


-5 92323 


-5 91695 


-5 92336 


-5.91213 


-5.92425 




(0.18352) 


(0.22155) 


(0.15302) 


(0.20352) 


(0.14326) 




1.60198 


1.60063 


1.60179 


1.60231 


1.61022 




(0.06385) 


(0.08602) 


(0.05271) 


(0.06851) 


(0.07535) 


Ho 


1.29726 


1.29592 


1.29800 


1.29429 


1.29766 




(0.09910) 


(0.12392) 


(0.08239) 


(0.10315) 


(0.09109) 




1.56533 


1.56363 


1.56472 


1.56333 


1.58932 




(0.09973) 


(0.12376) 


(0.08300) 


(0.11973) 


(0.12473) 


Ho 


1.03007 


1.02919 


1.03019 


1.03334 


1.03327 




(0.06328) 


(0.08554) 


(0.05223) 


(0.06728) 


(0.06998) 


0R 
Ho 


1.15416 


1.15318 


1.15416 


1.15399 


1.15456 




(0.06564) 


(0.08771) 


(0.05424) 


(0.07533) 


(0.06167) 


07 
Hi 


1.01944 


1.01857 


1.01991 


1.01744 


1.01934 




(0.06635) 


(0.08925) 


(0.05478) 


(0.06323) 


(0.06835) 


0s. 
Ho 


0.62225 


0.62171 


62038 


0.62111 


0.62448 




(0.06563) 


(0.08921) 


(0.05416) 


(0.07563) 


(0.07599) 


0Q 


56458 


56415 


56485 


56558 


56651 




(0.06018) 


(0.08321) 


(0.04957) 


(0.07066) 


(0.07538) 




0.35943 


0.35906 


0.36046 


0.35911 


0.36433 




(0.06714) 


(0.09243) 


(0.05532) 


(0.08332) 


(0.07755) 


Hll 


01035 


01034 


01035 


01035 


01037 




(0.00043) 


(0.00052) 


(0.00036) 


(0.00049) 


(0.00066) 


6i 


1.36408 


1.98699 


1.64216 


1.88568 


1.45548 




(1.22578) 


(1.22669) 


(1.22819) 


(1.22877) 


(1.22439) 


62 


0.01457 


0.01147 


0.01483 


0.01266 


0.01387 




(0.00361) 


(0.00362) 


(0.00362) 


(0.00374) 


(0.00369) 



estimates of /i, nevertheless, for the parameter (j), we note considerably differences between the maximum 
likelihood estimates and the corrected ones. 

8 Conclusions 

We defined a general beta regression models which allows a regression structure on the precision param- 
eter, and both the regression structures on the mean and on the precision parameters arc allowed to 
be nonlinear. Then, using the approximation theory developed by Cox and Snell (1968), we calculate 
O {n~^) bias for the MLEs for and 9. Our results generalize the formulae obtained by Ospina et al. 
(2006). We then defined bias-free estimator to order O (n~^'j, by using the expressions obtained through 
Cox and Snell's (1968) formulae and Firth's (1993) estimating equation. We also considered two schemes 
of bias correction based on bootstrap. We use simulation in a nonlinear beta regression model with 
nonlinear dispersion covariates to conclude the superiority of the corrective and parametric bootstrap 
methods of bias correction over the other methods presented here with regard to both bias reduction and 
mean squared error. In fact, the parametric bootstrap presented, in general, the smallest mean squared 
error, and the corrective method the smallest bias. Further, the corrective method has an advantage over 
the parametric bootstrap, which is the fact that the parametric bootstrap is computer intensive whereas 
the corrective method is not. The simulation also considered another study, to check if a nonlinear model 
which can be linearized should be estimated through a nonlinear mode;!, or through a linear model. The 
simulation showed that the nonlinear model should be preferred, showing that the general model we are 
presenting has the potential to be very useful to practioners. The paper is concluded with an empirical 
application to illustrate the usefulness of our results, and more specifBcally, the usefulness of considering 
a model with dispersion covariates. 
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Table 21: Estimated values of fi and 0. 



i 




MLE 

Ml 4'i 




Cox-Snell 




Firth 

Mi 4'i 




p-boot 

Mi 01 




nn- 


boot 

ft 


1 


0, 


.0999 


77.500 


0, 


.0999 


66 000 

\J\J . \J\J\J 


0, 


.0996 


92 800 


0, 


.0996 


66 900 


0. 


.0998 


53.600 


2 


0, 


.1865 


215.00 


0, 


.1866 


154.50 


0, 


.1863 


279.50 


0, 


.1863 


203 80 


0. 


.1863 


200.00 


3 


0, 


.3214 


596 30 


0, 


.3214 


311.10 


0, 


.3212 


798.10 


0, 


.3213 


577.10 


0. 


.3214 


658.30 


4 


0, 


.4737 


1471.7 


0, 


.4737 


432.20 


0, 


.4736 


1941.1 


0, 


.4737 


1371.2 


0. 


.4737 


1772.5 


5 


0, 


.0856 


93.700 


0, 


.0855 


77.900 


0, 


.0854 


114.40 


0, 


.0853 


82 800 


0. 


.0853 


69.400 


6 


0, 


.1421 


208 80 


0, 


.1419 


151.10 


0, 


.1418 


271.00 


0, 


.1418 


197.60 


0. 


.1419 


193.00 


7 


0, 


.2628 


61 3 QO 


0, 


.2626 


31 6 30 
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Appendix 

We give explicit expressions for the cumulants and their derivatives, both defined in pageO Further, we 
give the expressions for each quantity contained in equation ([S]) , some of them are also deduced to help 
the reader who might be interested in checking the results. 

We call attention for the fact that the expressions of the type contained in equation ^ from the 
Appendix of Ospina et al. (2006) are incorrect if we put a as cf). But this is a minor issue, and the 
formulae contained in the body of the text are correct. 

Consider initially the following notation for the derivatives, and product of the derivatives, of the 
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predictor with respect to the regression parameters: 



{rs)i 



-, iRS)i 



and so on. Now, consider the following quantities: 



, {rs,T)i 



•2i 



= V'((l-M»)</'i) + ^'(Mi<^i), 

= r{4>i) - rM)i^i - - Mi)<^i)(i - Mi)'- 



By using these quantities, the cumulants can be written as 

2 



KrS = - <l>i{lJ'iai - ^'((1 - S)i 



krs 



i=l 



dr]ii dr]2i 



dr}: 



•2i 



E, 2 j , ( dni\ djii d iiA 

9i { <PiCi -} — - ^ai- — } [r,s,u)i 

- (j>1ai {^^^ ■")» + (™' «)i + '')»}' 

" f d ■ \ ^ dd) 

^ \dr]uj dri2i 

+ X:M^'((l-Mi)<^.)-M.aO|^^(r.,t^)i, 
^ arm "'V2i 

nrsu = Em(l-M.)0O-M.aj|2|^(^) + |^ ^ | (r, 5, C/) 

+x:^.^.|^f^)'(-,5,f/)i 

+ f^,/.i{/((l-rt)<A0-Miai}^^(r-,5f/),, 
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Differentiating the second order cumulants with respect to the parameters, we have 

{{ru,s)i + {su,r)i}, 

= -i:{<^?[V'''((l-/^i)<^0-M.c.]+2</>,a,}(|^y^(r,s,C/),, 

41 = ER0i + 2^'((l-M.)<^,)-2M,aO^f^y(i?,5,«)„ 
^ artu \ari2ij 

_ V 6. (^) {{RU, S)i + {SU, R)i}, 

+ <l>iW{{l - - aiMO^^(r, S, U), 

■Hf dr]2i 
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We now define the diagonal matrices 



Ml = diag 



Ma 



drili 



' drjii d-ql^ 



aias f tec -<((!-, .)*.)} (I^l) * 



d'^m 



dVu dr}2t 



dni 



diag y — ^ I [2ai + c/j^ip ((1 - - ^iMiC^] [jj^ j 

d'^^i d(t>i 



d-ql^ d-q2i 



+ - ^JLi)(t)i) - flittii ^ 2 

diag f i J + 2^'i{l - - 2A*.a,]^ f ^ 



d^ij d'^(j)t 
drjii d-ql^ 



Ms = 



diag i 



dm 



' dri2i dr^l^ 



(20) 



(21) 



(22) 



(23) 
(24) 

(25) 



4>1ai ( d\x 



A^2 = diag ( ^[-0'((1 - - OiMilx^T^ ; ' 

2 dT^H «'72i / 



^3 - ^^^M^l^- 



(26) 
(27) 
(28) 



Let now, vfiji be the «th diagonal element of the matrix Mj, and riji be the ith diagonal element of 
the matrix iVj. Thus, consider also the following matrices 



L2 = ^ n2,Z,X^«X^,5„ 



i=l 



where 5^ is an n x 1 vector with a one in the zth position and zeros elsewhere. 
From expressions to (pS)) we get: 



2---" = ^rnij(r,s,u)j + ^ nij{(rs, u)^ - {ru,s)^ - (su,r)i}, 

^ n n 

'^rI ~ = ^m2i{R,s,u)i + ^n2i{su,R)i, 

4=1 i=l 
^ n n 

K^-S - ^'^'■Su = X! ™2»(?', S", it)i + ^ n2i{ru, S)i, 
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^ n n 

- -^KrsU = ^ "^3^(^', S, U)i - ^ n2j(rs, U)i, 
i=l 4=1 

-krsu = ^ rriiiiR, S, u)i - ^ n2i{RS, u)i, 



^KrStf = '^rn^i{r, S, U)i + '^n2i{SU,r)i, 

i=l i=l 
^ 71 n 



We are now in conditions to compute each term of expression l[8|: 



r ft 11 ^ 



1— 1 r s,u 
n 

+ X! '^1* X! K'"^K™{(rs, u)i - (ru, s)i - (su, r),} 

■i^l r,s,u 

n 

2— 1 r 

i—l r s,u 

n 
1=1 

= SlK^^X^MiPpp - S'^K'^'^X'^NiF, 

where X"'' is the matrix K^f^ if a = 1, . . . ,k and is K^^^ iia~k + l,...,k + h, further, if a = 
k + I, . . . , k + h, we use the abuse of notation 5a to mean Sa-h, also, the matrices Pf3f3, Fi and F were 
defined in pageO Similarly, we obtain, 

Y (41 - l^R^u] = S^K'^'Z^M2Pf3p + S^K^'Z^N^F, 

R,s,u ^ ^ 
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where, K"'^ is if^" if a = 1, . . . , fc and is K'^'^ if a ^ k + 1, . . . , k + h. Further, 

r,S,u ^ 1=1 r S.M 



+E E (^)'E'^'"'("^)» 

i=l S,u r 

n 

i—1 r 

n 

+ E "2'' E 

= <5fE^"''^'^'5'"^2.<5f(Xif'''^Z^)<5, 

+C E E i^-^^.'^.c^'''^^'^. 

i— 1 u 



In a similar fashion, we have 



where P^e was defined in page El Now, analogously we have, 



i '^flS ~ T^l^RSu > = O^K Z R'hPpe - d^K L2, 

E^aRsU \ ,AU) 1 1 xT i.^ae r, . xT T^^ae t 



where G was defined in page [HI Finally, we have 

R,S,U ^ ^ > 
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